Phase diagram of spin-1 bosons on one-dimensional lattices 
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Spinor Bose condensates loaded in optical lattices have a rich phase diagram characterized by 
different magnetic order. Here we apply the Density Matrix Renormalization Group to accurately 
determine the phase diagram for spin-1 bosons loaded on a one- dimensional lattice. The Mott lobes 
present an even or odd asymmetry associated to the boson filling. We show that for odd fillings 
the insulating phase is always in a dimerized state. The results obtained in this work are also 
relevant for the determination of the ground state phase diagram of the S — 1 Heisenberg model 
with biquadratic interaction. 
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The experimental realization of optical lattices Q has 
paved the way to study strongly correlated many-particle 
systems with cold atomic gases (see e.g. 0,13)- The main 
advantages with respect to condensed matter systems lie 
on the possibility of a precise knowledge of the under- 
lying microscopic models and an accurate and relatively 
easy control of the various couplings. Probably one of 
the most spectacular experiments in this respect is the 
observation Q of a Superfluid - Mott Insulator transi- 
tion previously predicted in [5( by a mapping onto the 
Bose-Hubbard model Q. 

More recently the use of far-off-resonance optical traps 
has opened the possibility to study spinor conden- 
sates pj. Spin effects are enhanced by the presence of 
strong interactions and small occupation number, thus 
resulting in a rich variety of phases with different mag- 
netic ordering. For spin-1 bosons it was predicted that 
the Mott insulating phases have nematic singlet [S| or 
dimerized || ground state depending on the mean oc- 
cupation and on the value of the spin exchange. Since 
the paper by Demler and Zhou Q several works have ad- 
dressed the properties of the phase diagram of spinor con- 
densates trapped in optical lattices (lrtllllll2Ul3ul4lll5l |. 
The increasing attention in spinor optical lattices has also 
revived the attention on open problems in the theory of 
quantum magnetism. The spinor Bose-Hubbard model, 
when the filling corresponds to one boson per site, can 
be mapped onto the S = 1 Heisenberg model with bi- 
quadratic interactions which exhibits a rich phase dia- 
gram including a lon g de bated nematic to dimer quantum 
phase transition HI lH GJ, H |H H |H H3 . 

Up to now the location of the phase boundary of 
the spinor Bose-Hubbard model has been determined by 
means of mean-field and strong coupling approaches. A 
quantitative calculation of the phase diagram is however 
still missing. This might be particularly important in 
one dimension where non-perturbative effects are more 
pronounced. This is the aim of this Letter. We deter- 
mine the location of the Mott lobes showing the even/odd 
asymmetry in the spinor case discussed in We then 
discuss the magnetic properties of the first lobe, conclud- 



ing that it is always in a dimerized phase. 

The effective Bose-Hubbard Hamiltonian, appropriate 
for S — 1 bosons, is given by 
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The operator a\ a creates a boson in the lowest Bloch 
band localized on site i and with spin component a 
along the quantization axis: hi = a\ a di tCr and Si = 

a i a\ ^Ta^'di^' are the total number of particles and 

the total spin on site i (T are the spin-1 operators). 
Atoms residing on the same lattice site have identical 
orbital wave function and their spin function must be 
symmetric. This constraint imposes that Si + rii must be 
even. The uniqueness of the completely symmetric state 
with fixed spin and number makes it possible to denote 
the single site states with \m, Si, S?). The coupling con- 
stants, which obey the constraint —1 < U2/U0 < 1/2, can 
be ex pre ssed in terms of the appropriate Wannier func- 
tions [l(J. Uo is set as the energy scale unit: Uq = 1. We 
discuss only the anti-ferromagnetic case (0 < U2 < 1/2). 

In the absence of spin dependent coupling a qualitative 
picture of the phase diagram can be drawn starting from 
the case of zero hopping (t — 0). The ground state is 
separated from any excited state by a finite energy gap. 
For finite hopping strength, the energy cost to add or re- 
move a particle AE± (excitation gap) is reduced and at 
a critical value ijr(/x) vanishes. This phase is named the 
Mott insulator. For large hopping amplitudes the ground 
state is a globally coherent superfluid phase. When U2 is 
different from zero, states with lowest spins, compatible 
with the constraint rii + Si = even, are favoured. This 
introduces an even/odd asymmetry of the lobes: the am- 
plitude of lobes with odd filling is reduced as compared 
with the lobes corresponding to even fillings 0- In the 
first lobe the extra energy required to have two particles 
on a site (instead of one) is 1 + 2J7 2 — /i, thus lower- 
ing the chemical potential value where the second lobe 
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starts. On the other hand, having no particles on a site 
gives no gain due to spin terms, accounting for the nearly 
unvaried bottom boundary of the lobe. 

In order to determine the phase diagram of Eq. Q we 
use the finite-size Density Matrix Renormalization Group 
(DMRG) with open boundary conditions |24| . The strat- 
egy of the DMRG is to construct a portion of the system 
(called the system block) and then recursively enlarge it, 
until the desired system size is reached. At every step 
the basis of the corresponding Hamiltonian is truncated, 
so that the size of the Hilbert space is kept manageable 
as the physical system grows. The truncation of the 
Hilbert space is performed by retaining the eigenstates 
corresponding to the m highest eigenvalues of the block's 
reduced density matrix. 

The DMRG has been employed, for the spinless case, 
in [2H l2t| . The presence of the spin degree of freedom 
makes the analysis considerably more difficult. In the 
numerical calculations the Hilbert space for the on-site 
part of the Hamiltonian is fixed by imposing a maximum 
occupation number n max . As the first lobe is character- 
ized by an insulating phase with n = 1 particle per site 
we choose n max = 3 in this case; the dimension of the 
Hilbert space per site becomes d = 20. We have checked, 
by increasing the value of n max , that this truncation of 
the Hilbert space is sufficient to compute the first lobe. 
In each DMRG iteration we keep up to m = 300 states 
in order to guarantee accurate results. The numerical 
calculations of the second lobe (n — 2 particles per site) 
have been performed with n max = 4 (which corresponds 
to d = 35). 

Phase Diagram - In the insulating phase the first ex- 
cited state is separated by the ground state by a Mott 
gap. In the limit of zero hopping the gap is determined 
by the extra energy AE± needed to place/remove a bo- 
son at a given site. The finite hopping renormalizes the 
gap which will vanish at a critical value. Then the sys- 
tem becomes supcrfhiid. This method has been emp loyed 
for the spinless case by Freericks and Monien |27|. and 



in |25|, 26] where it was combined with the DMRG. Here 
we use it for the spinor case. Three iterations of the 
DMRG procedure are performed, with projections on dif- 
ferent number sectors; the corresponding ground states 
give the desired energies E , E± = E + AE± . As target 
energies we used those obtained by the mapping of the 
Bose-Hubbard system into effective models as described 
in [ljj. We considered chains up to I = 128 sites for 
the first lobe, and L = 48 for the second lobe. The ex- 
trapolation procedure to extract the asymptotic values 
was obtained by means of linear fit in 1/L, as discussed 
in |26| . A comparison with a quadratic fit shows that 
0(1/ L 2 ) corrections are negligible on the scale of FigQ] 
The plot of the phase diagram in the t) plane for 
different values of the spin coupling U2 is shown in Fig^ 
The first lobe tends to reduce its size on increasing the 
spin coupling; in particular the upper critical chemical 
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FIG. 1: Phase diagram for the first two lobes of the ID 
Bose-Hubbard spin 1 model with nearest-neighbor interac- 
tions. The different panels correspond to different values of 
U2. The curves for U2 = coincide with the first two lobes 
for the spinless model computed in Refs.|25tl2q. 



potential at t = is M+(0) = 1 — 2U2, while the t* value 
of the hopping strength over which the system is always 
superfiuid is suppressed as U2 increases. On the other 
hand, the second lobe grows up when U2 increases. This 
even/odd effect, predicted in ||, is quantified in Fig^ 

Magnetic properties of the first Mott lobe - The first 
lobe of the spinor Bose lattice has a very interesting mag- 
netic structure. In presence of small hopping t boson tun- 
neling processes induce effective pairwise magnetic inter- 
actions between the spins, described by Hamiltonian |l0| : 



H cS = [cos 9 (Sj • Sj) + sinO (Sj • S^) 2 (2) 



with 



I 2t 2 

tan6»=- — k=- — Vl + tan 2 6> . (3) 



l-2U 2 



1 + U 2 



The absence of higher order terms, such as (S, ■ Sj) 3 , is 
due to the fact that the product of any three spin opera- 
tors can be expressed via lower order terms. In the case 
of anti- ferromagnetic interaction in Eq.Q, the parame- 
ter 9 varies in the interval 9 £ [— 3/47T, — tt/2[. Because of 
the form of the magnetic Hamiltonian, each bond tends 
to form a singlet-spin configuration, but singlet states on 
neighboring bonds are not allowed. There are two possi- 
ble ground states that may appear in this situation. A ne- 
matic state can be constructed by mixing states with to- 
tal spin S = and S — 2 on each bond. This construction 
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FIG. 2: Finite-size scaling of Dl for selected values of 9: 
circles (9 = — 0.657r), squares (— 0.7n), diamonds (— 0.727r), 
triangles up (— 0.737r), triangles down (— 0.7357r), triangles 
left (— 0.747t). In order to extrapolate the order parameter D, 
numerical data have been fitted with Dl = D+cL~ a (straight 
lines). DMRG simulations are performed with m ~ 140 for 
9 > -0.73tt, and m ~ 250 for 9 < -0.73tt. 



can be repeated on neighboring bonds, thereby preserv- 
ing translational invariance. This state breaks the spin- 
space rotational group 0(3), though time-reversal sym- 
metry is preserved. The expectation value of any spin op- 
erator vanishes ((Sf) =0, a = x,y, z), while some of the 
quadrupole operators have finite expectation values. The 
tensor Q ab = (S a S b ) - ^S ab is a traceless diagonal ma- 
trix, due to invariance under spin reflections. Since it has 
two identical eigenvalues «(5f ) 2 ) = ((Sf) 2 } ^ ((S'f) 2 }), 
it can be written as Q ab — Q (d a d b — |<5 ab ) using an or- 
der parameter (Q) = ((S?) 2 ) - <(Sf ) 2 ) = f((£f) 2 ) - 1 
and a unit vector d = ±z. However, since [Q,7i c tf] — 0, 
it is not possible to get Q ^ in finite-size systems, anal- 
ogously to what happens for the magnetization without 
external field. Therefore we characterized the range of 
nematic correlations in the ground state by coupling this 
operator to a fictitious "nematic field" : H\ = Hcs + XQ, 
and by evaluating the nematic susceptibility x ncm as a 
function of L: 



d 2 E (X) 



dA 2 
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(4) 



where E (\) is the ground energy of H\, Qo n is the 
matrix element between the ground and an excited state 
of H c ff (respectively with energy Eq and E 7 ). 

On the other hand a possibility to have SO(3) sym- 
metric solution stems from breaking translational invari- 
ance. Indeed a dimerized solution with singlets on every 
second bond satisfy these requirements. Dimerization 
could be described looking at the differences in expecta- 
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FIG. 3: Dimerization order parameter D near the ferromag- 
netic boundary: solid line shows a power law fit D ~ (9 — 9F) 1 
of numerical data with 7 ~ 6.15; dashed line shows an expo- 
nential law fit D ~ exp[-a/(0 - 9 F )~ 1/2 ] with a ~ 2.91. The 
linear fit is done over data for 9 < —0.7n, while the exponen- 
tial fit is for 9 < — 0.73-7T. DMRG calculations are performed 
with up to m ~ 300 states. Inset: extrapolated scaled gap 
A2-0 = (L — l)(i?2 — Eq) at the thermodynamic limit, close 
to 6 F . 



tion values of pair Hamiltonian Ti.^ on adjacent links 
(H e ff = J2(a) E3- The order parameter D reads 



D = 



L Gff 



H 



iff / 



(5) 



It has been proposed 0] that a narrow nematic re- 
gion exists between the ferromagnetic phase boundary 
{Of = — 37r/4, i.e. U2 — 0) and a critical angle 8c ~ 
— 0.7-7T (i.e. U2 ~ 10 -2 ), whereas a dimerized solution 
is favoured in the remaining anti-ferromagnetic region 
6c < < — 7r/2. This implies that the dimerization order 
parameter D should scale to zero in the whole nematic 
region. This possibility has been analyzed in Ref. 0] 
where it was suggested that D might go to exponen- 
tially near the ferromagnetic boundary, making it diffi- 
cult to detect the effective existence of the nematic phase. 
This interesting challenge has motivated numerical in- 
vestigations with different methods 
present new DMRG results which clarify the magnetic 
properties of the first Mott lobe (for sufficiently small 
hopping) and, consequently, of the biquadratic Heisen- 
berg chain. 

According to our numerical calculation there is no 
intermediate nematic phase, indeed we found a power 
law decay of the dimerization order parameter near 
Of = — Stt/A. The simulations of the bilinear-biquadratic 
model @ are less time and memory consuming than 
Bose-Hubbard ones, since the local Hilbert space has a 
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finite dimension d = 3. The number of block states kept 
during the renormalization procedure was chosen step 
by step in order to avoid artificial symmetry breaking. 
This careful treatment insures that there are no spurious 
sources of asymmetry like partially taking into account a 
probability multiplet. Here we considered up to m ~ 300 
states in order to obtain stable results. Raw numerical 
data are shown in FigEl where the finite-size dimeriza- 
tion parameter D(L) is plotted as a function of the chain 
length L (sec Eq. GD and |23). Finite-size scaling was 
used to extrapolate to the thermodynamic limit. After 
the extrapolation to the L — > oo limit, see Fig. we 
fitted the dimer order parameter with a power law 




where 7 - 6.1502 and do ~ 0.09177 ir (Fig. El solid line). 
We also tried to fit our data by an exponential law 

D = D e-WV^ 7 (7) 

as suggested in 0, with a ~ 2.911, D ~ 9.617; this fit 
seems to work for narrower regions (Fig. dashed line), 
however from our numerics we cannot exclude an expo- 
nential behavior of D in the critical region. The dimer- 
ized phase thus seems to survive up to the ferromagnetic 
phase boundary, independently from the chosen fitting 
form. This is also confirmed by the fact that the scaled 
gap between the ground state Eq and the lowest excited 
state E 2 (which is found to have total spin St = 2) seems 
not to vanish in the interesting region 9 > — 0.757T (see 
inset of Fig. [J). 

Moreover we analyzed the susceptibility of the chain to 
nematic ordering Xnem- The numerical data, presented 
in Fig. ^ show a power law behavior Xncm(L) oc L a as 
a function of the system size. The exponent a (shown 
in the inset) approaches the value a = 3a,s9^9p. 
This can also be confirmed by means of a perturbative 
calculation around the exact solution available at 9f', in- 
deed one obtains |Qo,7| 2 ~ L 2 and [E~ — Eq) ~ L _1 to 
be inserted in Eq. Q. The increase of the exponent for 
9 — > 9p indicates, as suggested in |2l|, that a tendency 
towards the nematic ordering is enhanced as the dimer 
order parameter goes to zero. 

Conclusions - In this Letter we analyzed, by means 
of a DMRG analysis, the phase diagram of the one- 
dimensional spinor boson condensate on an optical lat- 
tice. We determined quantitatively the shape of the first 
two Mott lobes, and the even/odd properties of the lobes. 
We furthermore discussed the magnetic properties of the 
first lobe. Our results indicate that the Mott insulator is 
always in a dimerized phase. 

This work was supported by IBM (2005 Faculty 
award), and by the European Community through grants 
RTNNANO, SQUBIT2. 
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FIG. 4: Nematic susceptibility Xnem as a function of the sys- 
tem size L. The various symbols refer to different values of 
6: circles (6 = — 0.657r), squares (— 0.771"), diamonds (— 0.737r), 
triangles up (— 0.74-71"), triangles down (— 0.7457r), triangles left 
(— 0.74757t). Straight lines show a power law fit Xnem = cL a 
of numerical data. Inset: exponent a as a function of 9. 
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